
function [c,ceq,Gc,Gceq]=eq_con_jac_thetas_NOlearn(pars,~,~,M,X,Z,R,sgi_nodes,sgi_weights,ep_sd,~)

c=[];

[N,T,K]=size(X);
KZ=size(Z,3);
Kth=size(R,2);

M=reshape(M,N*T,1);
R=repmat(R,T,1);

alpha=repmat(pars(1:N,1),T,1);
theta=pars(N+(1:2*Kth),1);
theta=reshape(theta,Kth,2);
beta=pars(N+2*Kth+(1:2*N),1);
v=pars(3*N+2*Kth+(1:N),1);
f=repmat(pars(4*N+2*Kth+(1:N),1),T,1);
% vs=pars(5*N+2*Kth+(1:N),1);
xi=pars(6*N+2*Kth+(1:K),1);
gamma=pars(6*N+2*Kth+K+(1:KZ),1);
kappa=pars((6*N+2*Kth+K+KZ+1):end,1);

ep_sd=repmat(ep_sd,T,1);

if nargout<=2
    
    ZZ=zeros(N,N);
    for kz=1:KZ
        ZZ=ZZ+Z(:,:,kz)*gamma(kz,1);
    end 
    Q=diag(v.*ep_sd(1:N,1))*exp(-ZZ)*diag(v.*ep_sd(1:N,1));
    P=kron(eye(2),tril(Q)+tril(Q,-1)');
    P=P\eye(2*N);
    
    XX=zeros(N,T);
    Pd=repmat(sqrt(diag(P\eye(2*N))),1,T);
    B=repmat(beta,1,T);
    for t=1:T
        XX(:,t)=permute(X(:,t,:),[1,3,2])*xi;
    end
    XX=reshape(XX,N*T,1);
    PdA=reshape(Pd(1:N,:),N*T,1);
    Pd=reshape(Pd(N+(1:N),:),N*T,1);
    BA=reshape(B(1:N,:),N*T,1);
    B=reshape(B(N+(1:N),:),N*T,1);

    UU0=exp(alpha+(R*theta(:,1)).*(PdA*sgi_nodes(:,1)'+BA+ep_sd*sgi_nodes(:,2)'));
    UU0=M.*(UU0./(1+UU0));
    UU1=exp(alpha+(R*theta(:,2)).*(Pd*sgi_nodes(:,1)'+B+ep_sd*sgi_nodes(:,2)')-f-XX-kappa);
    UU1=M.*(UU1./(1+UU1));

    ceq=(UU1-UU0)*sgi_weights;
    % ceq=[U1,U0];
    
else
    
    ZZ=zeros(N,N);
    for kz=1:KZ
        ZZ=ZZ+Z(:,:,kz)*gamma(kz);
    end 
    Q=diag(v.*ep_sd(1:N,1))*exp(-ZZ)*diag(v.*ep_sd(1:N,1));
    P=kron(eye(2),tril(Q)+tril(Q,-1)');
    P=P\eye(2*N);
    P=repmat(P,1,1,T);

    XX=zeros(N,T);
    Pd=repmat(sqrt(diag(P(:,:,1)\eye(2*N))),1,T);
    B=repmat(beta,1,T);
    for t=1:T
        XX(:,t)=permute(X(:,t,:),[1,3,2])*xi;
    end
    XX=reshape(XX,N*T,1);
    PdA=reshape(Pd(1:N,:),N*T,1);
    Pd=reshape(Pd(N+(1:N),:),N*T,1);
    BA=reshape(B(1:N,:),N*T,1);
    BD=reshape(B(N+(1:N),:),N*T,1);

    yA=PdA*sgi_nodes(:,1)'+BA+ep_sd*sgi_nodes(:,2)';
    yD=Pd*sgi_nodes(:,1)'+BD+ep_sd*sgi_nodes(:,2)';
    
    UU0=exp(alpha+(R*theta(:,1)).*yA);
    UU0=M.*(UU0./(1+UU0));
    UU1=exp(alpha+(R*theta(:,2)).*yD-f-XX-kappa);
    UU1=M.*(UU1./(1+UU1));

    ceq=(UU1-UU0)*sgi_weights;

    UU0=UU0.*(1-UU0);
    UU1=UU1.*(1-UU1);

    Gc=[];
    
    Gceq=zeros(N*T,length(pars));
    % alpha
    Gceq(:,1:N)=((UU1-UU0)*sgi_weights).*repmat(eye(N),T,1);
    % theta
    for k=1:Kth
        Gceq(:,N+k)=-(UU0.*R(:,k).*yA)*sgi_weights;
        Gceq(:,N+Kth+k)=(UU1.*R(:,k).*yD)*sgi_weights;
    end
    % betaA
    Gceq(:,N+2*Kth+(1:N))=(-(R*theta(:,1)).*UU0*sgi_weights).*repmat(eye(N),T,1);
    % betaD
    Gceq(:,N+2*Kth+N+(1:N))=((R*theta(:,2)).*UU1*sgi_weights).*repmat(eye(N),T,1);
    % v, gamma
    DS=zeros(2*N,N+1);
    for n=1:N
        DS(:,n+1)=DB_vg_NOlearn(Z,ZZ,v,ep_sd(1:N,1),n);
          % v
        Gceq(:,3*N+2*Kth+n)=((R*theta(:,2)).*UU1.*(repmat(DS(N+(1:N),n+1),T,1)*sgi_nodes(:,1)')-...
        (R*theta(:,1)).*UU0.*(repmat(DS(1:N,n+1),T,1)*sgi_nodes(:,1)'))*sgi_weights;
    end
    DS(:,1)=DB_vg_NOlearn(Z,ZZ,v,ep_sd(1:N,1),0);
    Gceq(:,6*N+2*Kth+K+(1:KZ))=((R*theta(:,2)).*UU1.*(repmat(DS(N+(1:N),1),T,1)*sgi_nodes(:,1)')-...
        (R*theta(:,1)).*UU0.*(repmat(DS(1:N,1),T,1)*sgi_nodes(:,1)'))*sgi_weights;
    % f
    A=-(UU1)*sgi_weights;
    Gceq(:,4*N+2*Kth+(1:N))=A.*repmat(eye(N),T,1);
    % kappa
    Gceq(:,(6*N+2*Kth+K+KZ+1):end)=diag(A);
    % xi
    for k=1:K
        Gceq(:,6*N+2*Kth+k)=-(UU1.*reshape(X(:,:,k),N*T,1))*sgi_weights;
    end
    
    Gceq=sparse(Gceq');
    
end


